/*
 Copyright (c) 2013 MeherTJ G. All rights reserved.
 License: LGPL for personnal study or free software.
 */

#pragma once

#include <memory>

#include "zen_point.h"

namespace Zen
{
	struct Mat2
	{
		union
		{
			float c[2][2]{};
			float v[4];
		} data;

		float& operator[](size_t i)
		{
			return data.v[i];
		}
		float const& operator[](size_t i) const
		{
			return data.v[i];
		}
	};

	struct Mat3
	{
		union
		{
			float c[3][3]{};
			float v[9];
		} data;

		float& operator[](size_t i)
		{
			return data.v[i];
		}
		float const& operator[](size_t i) const
		{
			return data.v[i];
		}
	};

	struct Mat4
	{
		union
		{
			float c[4][4]{};
			float v[16];
		} data;
		/*
		 // m30, m31, and m32 correspond to the translation values tx, ty, tz, respectively.
		 float m00, float m01, float m02, float m03,
		 float m10, float m11, float m12, float m13,
		 float m20, float m21, float m22, float m23,
		 float m30, float m31, float m32, float m33,
		 */

		float& operator[](size_t i)
		{
			return data.v[i];
		}
		float const& operator[](size_t i) const
		{
			return data.v[i];
		}
	};

	using ConstMat2 = Mat2 const;
	using ConstMat3 = Mat3 const;
	using ConstMat4 = Mat4 const;
} // namespace Zen

inline Zen::Mat3 operator+(Zen::ConstMat3& mL, Zen::ConstMat3& mR);
inline Zen::Mat3 operator-(Zen::ConstMat3& mL, Zen::ConstMat3& mR);
inline Zen::Mat3 operator*(Zen::ConstMat3& mL, Zen::ConstMat3& mR);
inline Zen::Mat3 operator/(Zen::ConstMat3& mL, Zen::ConstMat3& mR);
inline Zen::Mat4 operator+(Zen::ConstMat4& mL, Zen::ConstMat4& mR);
inline Zen::Mat4 operator-(Zen::ConstMat4& mL, Zen::ConstMat4& mR);
inline Zen::Mat4 operator*(Zen::ConstMat4& mL, Zen::ConstMat4& mR);
inline Zen::Mat4 operator/(Zen::ConstMat4& mL, Zen::ConstMat4& mR);
inline Zen::Mat4 operator&(Zen::ConstMat4& mL, Zen::ConstMat4& mR);

namespace Zen
{

	inline Mat2 Mat2Make(ConstMat3& matrix)
	{
		return { matrix[0], matrix[1], matrix[3], matrix[4] };
	}
	inline Mat2 Mat2Make(ConstMat4& matrix)
	{
		return { matrix[0], matrix[1], matrix[4], matrix[5] };
	}

	inline Mat3 Mat3Make(std::initializer_list<float> values)
	{
		Mat3 m3;
		auto dest = (float*)m3.data.v;
		for (auto i : values) *dest++ = i;
		return m3;
	}

	inline Mat3 Mat3Make(ConstMat4& matrix)
	{
		return { matrix[0], matrix[1], matrix[2], matrix[4], matrix[5], matrix[6], matrix[8], matrix[9], matrix[10] };
	}

	inline Mat3 Mat3MakeIdentity(float v = 1.f)
	{
		return Mat3{ v, 0, 0, 0, v, 0, 0, 0, v };
	}

	inline Mat3 Mat3MakeWithRows(ConstPoint3& row0, ConstPoint3& row1, ConstPoint3& row2)
	{
		return {
			row0[0], row1[0], row2[0], //
			row0[1], row1[1], row2[1], //
			row0[2], row1[2], row2[2], //
		};
	}

	inline Mat3 Mat3MakeWithColumns(ConstPoint3& col0, ConstPoint3& col1, ConstPoint3& col2)
	{
		return {
			col0[0], col0[1], col0[2], //
			col1[0], col1[1], col1[2], //
			col2[0], col2[1], col2[2], //
		};
	}

	inline Mat3 Mat3MakeScale(ConstPoint3& s)
	{
		return Mat3{ s.x, 0, 0, 0, s.y, 0, 0, 0, s.z };
	}

	inline Mat3 Mat3MakeRotation(float radians, ConstPoint3& axis)
	{
		Point3 v = Point3Normalize(axis);
		float cos = cosf(radians);
		float cosp = 1.0f - cos;
		float sin = sinf(radians);

		return {
			cos + cosp * v[0] * v[0],
			cosp * v[0] * v[1] + v[2] * sin,
			cosp * v[0] * v[2] - v[1] * sin,
			cosp * v[0] * v[1] - v[2] * sin,
			cos + cosp * v[1] * v[1],
			cosp * v[1] * v[2] + v[0] * sin,
			cosp * v[0] * v[2] + v[1] * sin,
			cosp * v[1] * v[2] - v[0] * sin,
			cos + cosp * v[2] * v[2],
		};
	}

	inline Mat3 Mat3MakeXRotation(float radians)
	{
		float cos = cosf(radians);
		float sin = sinf(radians);

		return { 1.0f, 0.0f, 0.0f, 0.0f, cos, sin, 0.0f, -sin, cos };
	}

	inline Mat3 Mat3MakeYRotation(float radians)
	{
		float cos = cosf(radians);
		float sin = sinf(radians);

		return { cos, 0.0f, -sin, 0.0f, 1.0f, 0.0f, sin, 0.0f, cos };
	}

	inline Mat3 Mat3MakeZRotation(float radians)
	{
		float cos = cosf(radians);
		float sin = sinf(radians);

		return { cos, sin, 0.0f, -sin, cos, 0.0f, 0.0f, 0.0f, 1.0f };
	}

	inline Mat2 Mat3GetMat2(ConstMat3& matrix)
	{
		Mat2 m = { matrix[0], matrix[1], matrix[3], matrix[4] };
		return m;
	}

	inline Point3 Mat3GetRow(ConstMat3& matrix, size_t row)
	{
		Point3 v = { matrix[row], matrix[3 + row], matrix[6 + row] };
		return v;
	}

	inline Point3 Mat3GetColumn(ConstMat3& matrix, size_t column)
	{
		Point3 v = { matrix[column * 3 + 0], matrix[column * 3 + 1], matrix[column * 3 + 2] };
		return v;
	}

	inline void Mat3SetRow(Mat3& matrix, size_t row, Point3 vector)
	{
		matrix[row] = vector[0];
		matrix[row + 3] = vector[1];
		matrix[row + 6] = vector[2];
	}

	inline void Mat3SetColumn(Mat3& matrix, size_t column, Point3 vector)
	{
		matrix[column * 3 + 0] = vector[0];
		matrix[column * 3 + 1] = vector[1];
		matrix[column * 3 + 2] = vector[2];
	}

	inline Mat3 Mat3Transpose(ConstMat3& o)
	{
		return { o[0], o[3], o[6], o[1], o[4], o[7], o[2], o[5], o[8] };
	}

	inline Mat3 Mat3Multiply(ConstMat3& mL, ConstMat3& mR)
	{
		Mat3 m;

		m[0] = mL[0] * mR[0] + mL[3] * mR[1] + mL[6] * mR[2];
		m[3] = mL[0] * mR[3] + mL[3] * mR[4] + mL[6] * mR[5];
		m[6] = mL[0] * mR[6] + mL[3] * mR[7] + mL[6] * mR[8];

		m[1] = mL[1] * mR[0] + mL[4] * mR[1] + mL[7] * mR[2];
		m[4] = mL[1] * mR[3] + mL[4] * mR[4] + mL[7] * mR[5];
		m[7] = mL[1] * mR[6] + mL[4] * mR[7] + mL[7] * mR[8];

		m[2] = mL[2] * mR[0] + mL[5] * mR[1] + mL[8] * mR[2];
		m[5] = mL[2] * mR[3] + mL[5] * mR[4] + mL[8] * mR[5];
		m[8] = mL[2] * mR[6] + mL[5] * mR[7] + mL[8] * mR[8];

		return m;
	}

	inline Mat3 Mat3Scale(ConstMat3& matrix, float x, float y, float z)
	{
		return {
			matrix[0] * x, matrix[1] * x, matrix[2] * x, //
			matrix[3] * y, matrix[4] * y, matrix[5] * y, //
			matrix[6] * z, matrix[7] * z, matrix[8] * z, //
		};
	}

	inline Mat3 Mat3Scale(ConstMat3& matrix, Point3 scalePointtor)
	{
		return {
			matrix[0] * scalePointtor[0], matrix[1] * scalePointtor[0], matrix[2] * scalePointtor[0], //
			matrix[3] * scalePointtor[1], matrix[4] * scalePointtor[1], matrix[5] * scalePointtor[1], //
			matrix[6] * scalePointtor[2], matrix[7] * scalePointtor[2], matrix[8] * scalePointtor[2], //
		};
	}

	inline Mat3 Mat3Scale(ConstMat3& matrix, Point4 scalePointtor)
	{
		return {
			matrix[0] * scalePointtor[0], matrix[1] * scalePointtor[0], matrix[2] * scalePointtor[0], //
			matrix[3] * scalePointtor[1], matrix[4] * scalePointtor[1], matrix[5] * scalePointtor[1], //
			matrix[6] * scalePointtor[2], matrix[7] * scalePointtor[2], matrix[8] * scalePointtor[2], //
		};
	}

	inline Mat3 Mat3Rotate(ConstMat3& matrix, float radians, ConstPoint3& axis)
	{
		Mat3 rm = Mat3MakeRotation(radians, axis);
		return Mat3Multiply(matrix, rm);
	}

	inline Mat3 Mat3RotateX(ConstMat3& matrix, float radians)
	{
		Mat3 rm = Mat3MakeXRotation(radians);
		return Mat3Multiply(matrix, rm);
	}

	inline Mat3 Mat3RotateY(ConstMat3& matrix, float radians)
	{
		Mat3 rm = Mat3MakeYRotation(radians);
		return Mat3Multiply(matrix, rm);
	}

	inline Mat3 Mat3RotateZ(ConstMat3& matrix, float radians)
	{
		Mat3 rm = Mat3MakeZRotation(radians);
		return Mat3Multiply(matrix, rm);
	}

	inline Point3 Mat3Multiply(ConstMat3& mL, Point3 vR)
	{
		return {
			mL[0] * vR[0] + mL[3] * vR[1] + mL[6] * vR[2],
			mL[1] * vR[0] + mL[4] * vR[1] + mL[7] * vR[2],
			mL[2] * vR[0] + mL[5] * vR[1] + mL[8] * vR[2],
		};
	}

	inline Mat4 Mat4Make(std::initializer_list<float> const& values)
	{
		Mat4 m;
		auto dest = (float*)m.data.v;
		for (auto i : values) *dest++ = i;
		return m;
	}
	inline Mat4 Mat4MakeIdentity(float v = 1.f)
	{
		return Mat4{ v, 0, 0, 0, 0, v, 0, 0, 0, 0, v, 0, 0, 0, 0, v };
	}

	inline Mat4 Mat4MakeWithRows(ConstPoint4& row0, ConstPoint4& row1, ConstPoint4& row2, ConstPoint4& row3)
	{
		return {
			row0[0], row1[0], row2[0], row3[0], //
			row0[1], row1[1], row2[1], row3[1], //
			row0[2], row1[2], row2[2], row3[2], //
			row0[3], row1[3], row2[3], row3[3], //
		};
	}

	inline Mat4 Mat4MakeWithColumns(ConstPoint4& col0, ConstPoint4& col1, ConstPoint4& col2, ConstPoint4& col3)
	{
		return {
			col0[0], col0[1], col0[2], col0[3], //
			col1[0], col1[1], col1[2], col1[3], //
			col2[0], col2[1], col2[2], col2[3], //
			col3[0], col3[1], col3[2], col3[3], //
		};
	}

	inline Mat4 Mat4MakeTranslation(ConstPoint3& t)
	{
		ConstMat4& m = {
			1, 0, 0, 0, //
			0, 1, 0, 0, //
			0, 0, 1, 0, //
			t.x, t.y, t.z, 1, //
		};
		return m;
	}

	inline Mat4 Mat4MakeScale(ConstPoint3& s)
	{
		return {
			s.x, 0, 0, 0, //
			0, s.y, 0, 0, //
			0, 0, s.z, 0, //
			0, 0, 0, 1, //
		};
	}

	/// Rotation begin
	inline Mat4 Mat4MakeRotation(float radians, ConstPoint3& r)
	{
		Point3 axis = Point3Normalize(r);
		float ct = cosf(radians);
		float st = sinf(radians);
		float ci = 1 - ct;
		float x = axis.x;
		float y = axis.y;
		float z = axis.z;

		return Mat4{
			ct + x * x * ci,
			y * x * ci + z * st,
			z * x * ci - y * st,
			0,
			x * y * ci - z * st,
			ct + y * y * ci,
			z * y * ci + x * st,
			0,
			x * z * ci + y * st,
			y * z * ci - x * st,
			ct + z * z * ci,
			0,
			0,
			0,
			0,
			1,
		};
	}

	inline Mat4 Mat4MakeXRotation(float radians)
	{
		float cos = cosf(radians);
		float sin = sinf(radians);

		return {
			1.0f, 0.0f, 0.0f, 0.0f, //
			0.0f, cos, sin, 0.0f, //
			0.0f, -sin, cos, 0.0f, //
			0.0f, 0.0f, 0.0f, 1.0f, //
		};
	}

	inline Mat4 Mat4MakeYRotation(float radians)
	{
		float cos = cosf(radians);
		float sin = sinf(radians);

		return {
			cos, 0.0f, -sin, 0.0f, //
			0.0f, 1.0f, 0.0f, 0.0f, //
			sin, 0.0f, cos, 0.0f, //
			0.0f, 0.0f, 0.0f, 1.0f, //
		};
	}

	inline Mat4 Mat4MakeZRotation(float radians)
	{
		float cos = cosf(radians);
		float sin = sinf(radians);

		return {
			cos, sin, 0.0f, 0.0f, //
			-sin, cos, 0.0f, 0.0f, //
			0.0f, 0.0f, 1.0f, 0.0f, //
			0.0f, 0.0f, 0.0f, 1.0f, //
		};
	}

	inline Mat4 Mat4MakePerspective(float fovyRadians, float aspect, float /*nearZ*/ n, float /*farZ*/ f)
	{
		float ys = 1 / tanf(fovyRadians * 0.5f);
		float xs = ys / aspect;
		float dz = n - f;
		float zs = (n + f) / dz;

		return {
			xs, 0, 0, 0, 0, ys, 0, 0, 0, 0, zs, -1, 0, 0, 2 * f * n / dz, 0 // 2 * f * n / dz - zs * n, +n,
		};

		/* move nearZ to zero.
		 return {
		 xs, 0, 0, 0,
		 0, ys, 0, 0,
		 0, 0, zs, -1,
		 0, 0, 2 * f * n / dz + 2 * f * n / dz - zs * n, +n,
		 };
		 */
	}

	inline Mat4 Mat4MakeFrustum(float left, float right, float bottom, float top, float nearZ, float farZ)
	{
		float ral = right + left;
		float rsl = right - left;
		float tsb = top - bottom;
		float tab = top + bottom;
		float fan = farZ + nearZ;
		float fsn = farZ - nearZ;

		return { 2.0f * nearZ / rsl, 0.0f, 0.0f, 0.0f, 0.0f, 2.0f * nearZ / tsb, 0.0f, 0.0f, ral / rsl, tab / tsb, -fan / fsn, -1.0f, 0.0f, 0.0f, (-2.0f * farZ * nearZ) / fsn, 0.0f };
	}

	inline Mat4 Mat4MakeOrtho(float left, float right, float bottom, float top, float nearZ, float farZ)
	{
		float ral = right + left;
		float rsl = right - left;
		float tab = top + bottom;
		float tsb = top - bottom;
		float fan = farZ + nearZ;
		float fsn = farZ - nearZ;

		return { 2.0f / rsl, 0.0f, 0.0f, 0.0f, 0.0f, 2.0f / tsb, 0.0f, 0.0f, 0.0f, 0.0f, -2.0f / fsn, 0.0f, -ral / rsl, -tab / tsb, -fan / fsn, 1.0f };
	}

	inline Mat4 Mat4MakeLookAt(ConstPoint3& ev, ConstPoint3& cv, ConstPoint3& uv)
	{
		Point3 n = Point3Normalize(ev - cv);
		Point3 u = Point3Normalize(Point3CrossProduct(uv, n));
		Point3 v = Point3CrossProduct(n, u);

		return { u[0], v[0], n[0], 0.0f, u[1], v[1], n[1], 0.0f, u[2], v[2], n[2], 0.0f, Point3DotProduct(-(u), ev), Point3DotProduct(-(v), ev), Point3DotProduct(-(n), ev), 1.0f };
	}

	inline Point4 Mat4GetRow(ConstMat4& matrix, size_t row)
	{
		Point4 v = { matrix[row], matrix[4 + row], matrix[8 + row], matrix[12 + row] };
		return v;
	}

	inline Point4 Mat4GetColumn(ConstMat4& matrix, size_t column)
	{
		Point4 v = { matrix[column * 4 + 0], matrix[column * 4 + 1], matrix[column * 4 + 2], matrix[column * 4 + 3] };
		return v;
	}

	inline void Mat4SetRow(Mat4& matrix, size_t row, ConstPoint4& v)
	{
		matrix[row] = v[0];
		matrix[row + 4] = v[1];
		matrix[row + 8] = v[2];
		matrix[row + 12] = v[3];
	}

	inline void Mat4SetColumn(Mat4& matrix, size_t column, ConstPoint4& v)
	{
		matrix[column * 4 + 0] = v[0];
		matrix[column * 4 + 1] = v[1];
		matrix[column * 4 + 2] = v[2];
		matrix[column * 4 + 3] = v[3];
	}

	/// multiply
	inline Point4 Mat4Multiply(ConstMat4& mL, ConstPoint4& vR)
	{
		return { mL[0] * vR[0] + mL[4] * vR[1] + mL[8] * vR[2] + mL[12] * vR[3], mL[1] * vR[0] + mL[5] * vR[1] + mL[9] * vR[2] + mL[13] * vR[3],
			mL[2] * vR[0] + mL[6] * vR[1] + mL[10] * vR[2] + mL[14] * vR[3], mL[3] * vR[0] + mL[7] * vR[1] + mL[11] * vR[2] + mL[15] * vR[3] };
	}

	inline Point3 Mat4Multiply(ConstMat4& mL, ConstPoint3& vR)
	{
		Point4 v4 = Mat4Multiply(mL, Point4Make(vR, 0.f));
		return Point3Make(v4[0], v4[1], v4[2]);
	}

	inline Point3 Mat4MultiplyAndProject(ConstMat4& mL, ConstPoint3& vR)
	{
		Point4 v4 = Mat4Multiply(mL, Point4Make(vR, 1.0f));
		return Point3Make(v4[0], v4[1], v4[2]) * (1.0f / v4[3]);
	}

	inline Mat4 Mat4Multiply(ConstMat4& mL, ConstMat4& mR)
	{
		Mat4 m;

		m[0] = mL[0] * mR[0] + mL[4] * mR[1] + mL[8] * mR[2] + mL[12] * mR[3];
		m[4] = mL[0] * mR[4] + mL[4] * mR[5] + mL[8] * mR[6] + mL[12] * mR[7];
		m[8] = mL[0] * mR[8] + mL[4] * mR[9] + mL[8] * mR[10] + mL[12] * mR[11];
		m[12] = mL[0] * mR[12] + mL[4] * mR[13] + mL[8] * mR[14] + mL[12] * mR[15];

		m[1] = mL[1] * mR[0] + mL[5] * mR[1] + mL[9] * mR[2] + mL[13] * mR[3];
		m[5] = mL[1] * mR[4] + mL[5] * mR[5] + mL[9] * mR[6] + mL[13] * mR[7];
		m[9] = mL[1] * mR[8] + mL[5] * mR[9] + mL[9] * mR[10] + mL[13] * mR[11];
		m[13] = mL[1] * mR[12] + mL[5] * mR[13] + mL[9] * mR[14] + mL[13] * mR[15];

		m[2] = mL[2] * mR[0] + mL[6] * mR[1] + mL[10] * mR[2] + mL[14] * mR[3];
		m[6] = mL[2] * mR[4] + mL[6] * mR[5] + mL[10] * mR[6] + mL[14] * mR[7];
		m[10] = mL[2] * mR[8] + mL[6] * mR[9] + mL[10] * mR[10] + mL[14] * mR[11];
		m[14] = mL[2] * mR[12] + mL[6] * mR[13] + mL[10] * mR[14] + mL[14] * mR[15];

		m[3] = mL[3] * mR[0] + mL[7] * mR[1] + mL[11] * mR[2] + mL[15] * mR[3];
		m[7] = mL[3] * mR[4] + mL[7] * mR[5] + mL[11] * mR[6] + mL[15] * mR[7];
		m[11] = mL[3] * mR[8] + mL[7] * mR[9] + mL[11] * mR[10] + mL[15] * mR[11];
		m[15] = mL[3] * mR[12] + mL[7] * mR[13] + mL[11] * mR[14] + mL[15] * mR[15];

		return m;
	}
	/// multiply end

	inline Mat4 Mat4Transpose(ConstMat4& m)
	{
		return { m[0], m[4], m[8], m[12], m[1], m[5], m[9], m[13], m[2], m[6], m[10], m[14], m[3], m[7], m[11], m[15] };
	}

	inline Mat4 Mat4Translate(ConstMat4& matrix, ConstPoint3& t)
	{
		return { matrix[0], matrix[1], matrix[2], matrix[3], matrix[4], matrix[5], matrix[6], matrix[7], matrix[8], matrix[9], matrix[10], matrix[11],
			matrix[0] * t[0] + matrix[4] * t[1] + matrix[8] * t[2] + matrix[12], matrix[1] * t[0] + matrix[5] * t[1] + matrix[9] * t[2] + matrix[13],
			matrix[2] * t[0] + matrix[6] * t[1] + matrix[10] * t[2] + matrix[14], matrix[3] * t[0] + matrix[7] * t[1] + matrix[11] * t[2] + matrix[15] };
	}

	inline Mat4 Mat4Translate(ConstMat4& matrix, ConstPoint4& t)
	{
		return {
			matrix[0],
			matrix[1],
			matrix[2],
			matrix[3],
			matrix[4],
			matrix[5],
			matrix[6],
			matrix[7],
			matrix[8],
			matrix[9],
			matrix[10],
			matrix[11],
			matrix[0] * t[0] + matrix[4] * t[1] + matrix[8] * t[2] + matrix[12],
			matrix[1] * t[0] + matrix[5] * t[1] + matrix[9] * t[2] + matrix[13],
			matrix[2] * t[0] + matrix[6] * t[1] + matrix[10] * t[2] + matrix[14],
			matrix[3] * t[0] + matrix[7] * t[1] + matrix[11] * t[2] + matrix[15],
		};
	}

	inline Mat4 Mat4Scale(ConstMat4& matrix, ConstPoint3& scale)
	{
		return {
			matrix[0] * scale[0],
			matrix[1] * scale[0],
			matrix[2] * scale[0],
			matrix[3] * scale[0],
			matrix[4] * scale[1],
			matrix[5] * scale[1],
			matrix[6] * scale[1],
			matrix[7] * scale[1],
			matrix[8] * scale[2],
			matrix[9] * scale[2],
			matrix[10] * scale[2],
			matrix[11] * scale[2],
			matrix[12],
			matrix[13],
			matrix[14],
			matrix[15],
		};
	}

	inline Mat4 Mat4Rotate(ConstMat4& matrix, float radians, Point3 axis)
	{
		Mat4 rm = Mat4MakeRotation(radians, axis);
		return Mat4Multiply(matrix, rm);
	}

	inline Mat4 Mat4RotateX(ConstMat4& matrix, float radians)
	{
		Mat4 rm = Mat4MakeXRotation(radians);
		return Mat4Multiply(matrix, rm);
	}

	inline Mat4 Mat4RotateY(ConstMat4& matrix, float radians)
	{
		Mat4 rm = Mat4MakeYRotation(radians);
		return Mat4Multiply(matrix, rm);
	}

	inline Mat4 Mat4RotateZ(ConstMat4& matrix, float radians)
	{
		Mat4 rm = Mat4MakeZRotation(radians);
		return Mat4Multiply(matrix, rm);
	}

	inline Point4 QuatNormalize(Point4 quat)
	{
		auto dis = Point4Length(quat);
		return quat / dis;
	}
	inline Mat3 Mat3MakeWithQuat(Point4 const& q)
	{
		auto qua = QuatNormalize(q);

		float x = qua[0];
		float y = qua[1];
		float z = qua[2];
		float w = qua[3];

		float _2x = x + x;
		float _2y = y + y;
		float _2z = z + z;
		float _2w = w + w;

		return {
			1.0f - _2y * y - _2z * z,
			_2x * y + _2w * z,
			_2x * z - _2w * y,
			_2x * y - _2w * z,
			1.0f - _2x * x - _2z * z,
			_2y * z + _2w * x,
			_2x * z + _2w * y,
			_2y * z - _2w * x,
			1.0f - _2x * x - _2y * y,
		};
	}

	inline Mat4 Mat4MakeWithQuat(Point4 qua)
	{
		qua = QuatNormalize(qua);

		float x = qua[0];
		float y = qua[1];
		float z = qua[2];
		float w = qua[3];

		float _2x = x + x;
		float _2y = y + y;
		float _2z = z + z;
		float _2w = w + w;

		return {
			1.0f - _2y * y - _2z * z,
			_2x * y + _2w * z,
			_2x * z - _2w * y,
			0.0f,
			_2x * y - _2w * z,
			1.0f - _2x * x - _2z * z,
			_2y * z + _2w * x,
			0.0f,
			_2x * z + _2w * y,
			_2y * z - _2w * x,
			1.0f - _2x * x - _2y * y,
			0.0f,
			0.0f,
			0.0f,
			0.0f,
			1.0f,
		};
	}

} // namespace Zen

inline Zen::Mat3 operator+(Zen::ConstMat3& mL, Zen::ConstMat3& mR)
{
	return {
		mL[0] + mR[0],
		mL[1] + mR[1],
		mL[2] + mR[2],
		mL[3] + mR[3],
		mL[4] + mR[4],
		mL[5] + mR[5],
		mL[6] + mR[6],
		mL[7] + mR[7],
		mL[8] + mR[8],
	};
}

inline Zen::Mat3 operator-(Zen::ConstMat3& mL, Zen::ConstMat3& mR)
{
	return {
		mL[0] - mR[0],
		mL[1] - mR[1],
		mL[2] - mR[2],
		mL[3] - mR[3],
		mL[4] - mR[4],
		mL[5] - mR[5],
		mL[6] - mR[6],
		mL[7] - mR[7],
		mL[8] - mR[8],
	};
}

inline Zen::Mat3 operator*(Zen::ConstMat3& mL, Zen::ConstMat3& mR)
{
	return {
		mL[0] * mR[0],
		mL[1] * mR[1],
		mL[2] * mR[2],
		mL[3] * mR[3],
		mL[4] * mR[4],
		mL[5] * mR[5],
		mL[6] * mR[6],
		mL[7] * mR[7],
		mL[8] * mR[8],
	};
}

inline Zen::Mat3 operator/(Zen::ConstMat3& mL, Zen::ConstMat3& mR)
{
	return {
		mL[0] / mR[0],
		mL[1] / mR[1],
		mL[2] / mR[2],
		mL[3] / mR[3],
		mL[4] / mR[4],
		mL[5] / mR[5],
		mL[6] / mR[6],
		mL[7] / mR[7],
		mL[8] / mR[8],
	};
}

inline Zen::Mat4 operator+(Zen::ConstMat4& mL, Zen::ConstMat4& mR)
{
	return {
		mL[0] + mR[0],
		mL[1] + mR[1],
		mL[2] + mR[2],
		mL[3] + mR[3],
		mL[4] + mR[4],
		mL[5] + mR[5],
		mL[6] + mR[6],
		mL[7] + mR[7],
		mL[8] + mR[8],
		mL[9] + mR[9],
		mL[10] + mR[10],
		mL[11] + mR[11],
		mL[12] + mR[12],
		mL[13] + mR[13],
		mL[14] + mR[14],
		mL[15] + mR[15],
	};
}

inline Zen::Mat4 operator-(Zen::ConstMat4& mL, Zen::ConstMat4& mR)
{
	return {
		mL[0] - mR[0],
		mL[1] - mR[1],
		mL[2] - mR[2],
		mL[3] - mR[3],
		mL[4] - mR[4],
		mL[5] - mR[5],
		mL[6] - mR[6],
		mL[7] - mR[7],
		mL[8] - mR[8],
		mL[9] - mR[9],
		mL[10] - mR[10],
		mL[11] - mR[11],
		mL[12] - mR[12],
		mL[13] - mR[13],
		mL[14] - mR[14],
		mL[15] - mR[15],
	};
}

inline Zen::Mat4 operator*(Zen::ConstMat4& mL, Zen::ConstMat4& mR)
{
	return {
		mL[0] * mR[0],
		mL[1] * mR[1],
		mL[2] * mR[2],
		mL[3] * mR[3],

		mL[4] * mR[4],
		mL[5] * mR[5],
		mL[6] * mR[6],
		mL[7] * mR[7],

		mL[8] * mR[8],
		mL[9] * mR[9],
		mL[10] * mR[10],
		mL[11] * mR[11],

		mL[12] * mR[12],
		mL[13] * mR[13],
		mL[14] * mR[14],
		mL[15] * mR[15],
	};
}

inline Zen::Mat4 operator/(Zen::ConstMat4& mL, Zen::ConstMat4& mR)
{
	return {
		mL[0] / mR[0],
		mL[1] / mR[1],
		mL[2] / mR[2],
		mL[3] / mR[3],
		mL[4] / mR[4],
		mL[5] / mR[5],
		mL[6] / mR[6],
		mL[7] / mR[7],
		mL[8] / mR[8],
		mL[9] / mR[9],
		mL[10] / mR[10],
		mL[11] / mR[11],
		mL[12] / mR[12],
		mL[13] / mR[13],
		mL[14] / mR[14],
		mL[15] / mR[15],
	};
}

inline Zen::Mat4 operator&(Zen::ConstMat4& mL, Zen::ConstMat4& mR)
{
	return Zen::Mat4Multiply(mL, mR);
}

inline Zen::Point4 operator&(Zen::ConstMat4& mL, Zen::ConstPoint4& vR)
{
	return Zen::Mat4Multiply(mL, vR);
}

template <typename CharT>
inline std::basic_ostream<CharT>& operator<<(std::basic_ostream<CharT>& o, Zen::ConstMat3& m)
{
	for (int i = 0; i < 3; ++i)
	{
		o << m.data.c[i][0];
		for (int j = 1; j < 3; ++j)
		{
			o << ',' << m.data.c[i][j];
		}
		o << '\n';
	}
	return o;
}

template <typename CharT>
inline std::basic_ostream<CharT>& operator<<(std::basic_ostream<CharT>& o, Zen::ConstMat4& m)
{
	using namespace std;
	for (int i = 0; i < 4; ++i)
	{
		o << std::to_string(m.data.c[i][0]);
		for (int j = 1; j < 4; ++j)
		{
			o << ',' << std::to_string(m.data.c[i][j]);
		}
		o << '\n';
	}
	return o;
}
